Statistical Rethinking (Code)

Code from Statistical Rethinking modified by R Pruim is shown below. Differences to the oringal include:

• a preference for putting data into containers (data frames, mostly), rather than working with lose vectors.
• use of lattice and/or ggplot2 rather than base graphics
• use of tidyverse for data transformation
• better (in my opinion) naming conventions

Intro: Waffle Houses and Divorce

The following model suggests that divorce rates are higher in states with more Waffle Houses (per capita).

data(WaffleDivorce)
WaffleDivorce <- 
  WaffleDivorce %>% 
  mutate(
    WaffleHouses.pm = WaffleHouses / Population
  )
m5.0 <- map(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu <- a + bW * WaffleHouses.pm,
    a ~ dnorm(10, 10),
    bW ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = WaffleDivorce
)
precis(m5.0, digits = 3)

##        Mean StdDev  5.5% 94.5%
## a     9.321  0.272 8.887 9.755
## bW    0.074  0.027 0.032 0.117
## sigma 1.677  0.168 1.409 1.945

Here is a plot.

require(statisticalModeling)   # need the version from rpruim on github
coef(m5.0)

##          a         bW      sigma 
## 9.32060116 0.07433683 1.67730153

gf_point(Divorce ~ WaffleHouses.pm, data = WaffleDivorce, color = rangi2) %>%
  gf_coefline(coef = coef(m5.0))

## Warning in gf_coefline(., coef = coef(m5.0)): Ignoring all but first two
## values of coef.

R code 5.1

# load data
library(rethinking)
data(WaffleDivorce)

# standardize predictor
WaffleDivorce <- 
  WaffleDivorce %>% 
  mutate(MedianAgeMarriage.s = zscore(MedianAgeMarriage))

# fit model
m5.1 <- map(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu <- a + bA * MedianAgeMarriage.s,
    a ~ dnorm(10, 10),
    bA ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = WaffleDivorce
)

R code 5.2

# compute percentile interval of mean
m5.1.pred <-
  data.frame(
    MedianAgeMarriage.s = seq(from = -3, to = 3, by = 0.25)
  )
m5.1.link <- link(m5.1, data = m5.1.pred)

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

m5.1.pred <- 
  m5.1.pred %>% 
  mutate(
    mu.PIlo = apply(m5.1.link, 2, PI)[1,],
    mu.PIhi = apply(m5.1.link, 2, PI)[2,]
  )

# plot it all
gf_point(Divorce ~ MedianAgeMarriage.s + label:Loc, color = rangi2, 
         data = WaffleDivorce) %>% 
  gf_ribbon(mu.PIlo + mu.PIhi ~ MedianAgeMarriage.s, data = m5.1.pred,  
            alpha = 0.2) %>%
  gf_coefline(coef = coef(m5.1), col = "red", alpha = 0.5) %>%
plotly::ggplotly()

## Warning: Ignoring unknown aesthetics: label

## Warning in gf_coefline(., coef = coef(m5.1), col = "red", alpha = 0.5):
## Ignoring all but first two values of coef.

R code 5.3

WaffleDivorce <-
  WaffleDivorce %>% 
  mutate(Marriage.s = zscore(Marriage))

m5.2 <- map(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu <- a + bR * Marriage.s,
    a ~ dnorm(10, 10),
    bR ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = WaffleDivorce
)

m5.2.pred <-
  data.frame(
    Marriage.s = seq(from = -3, to = 3, by = 0.25)
  )
m5.2.link <- link(m5.2, data = m5.2.pred)

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

m5.2.pred <- 
  m5.2.pred %>% 
  mutate(
    mu.PIlo = apply(m5.2.link, 2, PI)[1,],
    mu.PIhi = apply(m5.2.link, 2, PI)[2,]
  )

# plot it all
gf_point(Divorce ~ Marriage.s, data = WaffleDivorce, col = rangi2) %>%
  gf_coefline(coef = coef(m5.2), col = "red", alpha = 0.5) %>%
  gf_ribbon(mu.PIlo + mu.PIhi ~ Marriage.s, data = m5.2.pred, alpha = 0.1)

## Warning in gf_coefline(., coef = coef(m5.2), col = "red", alpha = 0.5):
## Ignoring all but first two values of coef.

plotly::ggplotly()

R code 5.4

m5.3 <- map(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu <- a + bR * Marriage.s + bA * MedianAgeMarriage.s,
    a ~ dnorm(10, 10),
    bR ~ dnorm(0, 1),
    bA ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = WaffleDivorce
)
precis(m5.3)

##        Mean StdDev  5.5% 94.5%
## a      9.69   0.20  9.36 10.01
## bR    -0.13   0.28 -0.58  0.31
## bA    -1.13   0.28 -1.58 -0.69
## sigma  1.44   0.14  1.21  1.67

R code 5.5

plot(precis(m5.3))

R code 5.6

m5.4 <- map(
  alist(
    Marriage.s ~ dnorm(mu, sigma),
    mu <- a + b * MedianAgeMarriage.s,
    a ~ dnorm(0, 10),
    b ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = WaffleDivorce
)

R code 5.7 – computing residuals (the hard way)

# for each state, compute expected value at MAP and residual
WaffleDivorce <-
  WaffleDivorce %>% 
  mutate(
    mu.m5.4 = coef(m5.4)['a'] + coef(m5.4)['b'] * MedianAgeMarriage.s,
    resid.m5.4 = Marriage.s - mu.m5.4
  )

R code 5.7a – computing residuals (the “easy” way)

link() takes care of computing the model value for each posterior sample so we don’t have to access the coefficients directly. As models become more complicated, this will be easier. Also, it provides us with more information, since we have a distribution of model responses for each case.

We are using the defaults here, but link() takes two additional arguments in addition to the model.

• n: number of posterior samples to use. (The default is 1000.) Increase this to reduce sampling variability (at the cost of time).
• data: optional data set to use in place of original data from map(). Use this to make predictions using any predictor values you like.

# for each state, compute expected value at MAP and residual
m5.4.link <- link(m5.4) 

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

WaffleDivorce <-
  WaffleDivorce %>%
  mutate(
    mu = apply(m5.4.link, 2, mean),
    resid = Marriage.s - mu
  )
WaffleDivorce %>% select(Loc, mu, resid) %>% head(3)

##   Loc        mu       resid
## 1  AL 0.4289818 -0.40633776
## 2  AK 0.4858405  1.06396115
## 3  AZ 0.1446885 -0.09571417

R code 5.8

gf_segment(mu.m5.4 + Marriage.s ~ MedianAgeMarriage.s + MedianAgeMarriage.s, 
           data = WaffleDivorce, color = "red", alpha = 0.5) %>%
  gf_coefline(coef = coef(m5.4)) %>%
  gf_point(Marriage.s ~ MedianAgeMarriage.s, col = rangi2, data = WaffleDivorce) 

## Warning in gf_coefline(., coef = coef(m5.4)): Ignoring all but first two
## values of coef.

plotly::ggplotly()

R code 5.9

# prepare new counterfactual data
m5.3.pred <- 
  data_frame(
    Marriage.s = seq(from = -3, to = 3, by = 0.25),
    MedianAgeMarriage.s = mean(WaffleDivorce$MedianAgeMarriage.s)
  )

# compute counterfactual mean divorce (mu)
mu <- link(m5.3, data = m5.3.pred)

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

# simulate counterfactual divorce outcomes
R.sim <- sim(m5.3, data = m5.3.pred, n = 1e4)

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

m5.3.pred <-
  m5.3.pred %>% 
  mutate(
     mu.mean = apply(mu, 2, mean),
     mu.PIlo =  apply(mu, 2, PI)[1,],
     mu.PIhi =  apply(mu, 2, PI)[2,],
     R.PIlo = apply(R.sim, 2, PI)[1,],
     R.PIhi = apply(R.sim, 2, PI)[2,]
  )

gf_line(mu.mean ~ Marriage.s, data = m5.3.pred, color = "gray50") %>%
  gf_ribbon(mu.PIlo + mu.PIhi ~ Marriage.s, fill = "gray50", alpha = 0.2) %>%
  gf_ribbon(R.PIlo + R.PIhi ~ Marriage.s, fill = "gray50", alpha = 0.2) +
  labs( y = "Divorce rate", caption = "MedianAgeMarriage.s = 0")

R code 5.10

m5.3.pred2 <-
  data_frame(
    MedianAgeMarriage.s =  seq(from = -3, to = 3, by = 0.25),
    Marriage.s = mean(WaffleDivorce$Marriage.s)
    )

mu <- link(m5.3, data = m5.3.pred2, n = 1e4)

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

A.sim <- sim(m5.3, data = m5.3.pred2, n = 1e4)

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

m5.3.pred2 <-
  m5.3.pred2 %>%
  mutate(
    mu.mean = apply(mu, 2, mean),
    mu.PIlo = apply(mu, 2, PI)[1,],
    mu.PIhi = apply(mu, 2, PI)[2,],
    A.PIlo = apply(A.sim, 2, PI)[1,],
    A.PIhi = apply(A.sim, 2, PI)[2,]
  )

gf_line(mu.mean ~ MedianAgeMarriage.s, data = m5.3.pred2, color = "red") %>%
  gf_ribbon(mu.PIlo + mu.PIhi ~ MedianAgeMarriage.s, alpha = 0.1) %>%
  gf_ribbon(A.PIlo + A.PIhi ~ MedianAgeMarriage.s, alpha = 0.1) +
  labs( y = "Divorce rate", caption = "Marriage.s = 0")

R code 5.11

# call link without specifying new data
# so it uses original data
divorce.mu  <- link(m5.3)

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

# simulate observations
# again no new data, so uses original data
divorce.sim <- sim(m5.3, n = 1e4)

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

# summarize samples across cases
m5.3.pred3 <-
  WaffleDivorce %>% 
  mutate(
    mu.mean = apply(divorce.mu, 2, mean),
    mu.PIlo = apply(divorce.mu, 2, PI)[1,],
    mu.PIhi = apply(divorce.mu, 2, PI)[2,],
    divorce.PIlo = apply(divorce.sim, 2, PI)[1,],
    divorce.PIhi = apply(divorce.sim, 2, PI)[2,]
  )

R code 5.12

gf_pointrange(mu.mean + divorce.PIlo + divorce.PIhi ~ Divorce, data = m5.3.pred3,
              color = rangi2) %>%
  gf_abline(slope = 1, intercept = 0) +
  labs(x =  "Observed divorce", y = "Predicted divorce")

R code 5.13

# make most recent ggplot interactive
# adding text to the geom_point() above makes it available on hover here.
plotly::ggplotly()

R code 5.14

# compute residuals
m5.3.pred3 <-
  m5.3.pred3 %>% 
  mutate(divorce.resid = Divorce - mu.mean,
         state = reorder(Loc, divorce.resid))

gf_linerange((Divorce - divorce.PIlo) + (Divorce - divorce.PIhi) ~ Loc, 
             data = m5.3.pred3, alpha = 0.2, size = 0.9) %>%
  gf_linerange((Divorce - mu.PIlo) + (Divorce - mu.PIhi) ~ Loc, 
             data = m5.3.pred3, alpha = 0.2, size = 1.4) %>%
  gf_point(divorce.resid ~ Loc, data = m5.3.pred3) + 
  coord_flip() 

WaffleDivorce <-
  WaffleDivorce %>%
  mutate(
    WaffleHouses.pc = WaffleHouses / Population * 1e6,
    divorce.resid = Divorce - apply(link(m5.3), 2, mean) 
  )

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

gf_point(divorce.resid ~ WaffleHouses.pc + label:Loc, 
         data = WaffleDivorce) %>%
  gf_smooth(divorce.resid ~ WaffleHouses.pc, method = "lm", size = 0.5) %>%
  plotly::ggplotly()

## Warning: Ignoring unknown aesthetics: label

### R code 5.15

n <- 100
Sim5.15Data <-
  data_frame( 
    x_real = rnorm(n),         # x_real as Gaussian with mean 0 and stddev 1
    x_spur = rnorm(n, x_real), # x_spur as Gaussian with mean = x_real
    y =  rnorm(n, x_real)      # y as Gaussian with mean=x_real
  )

m_both <-
  map(
    alist(
      y ~ dnorm(mu, sigma),
      mu <- a + b_real * x_real + b_spur * x_spur,
      a ~ dnorm(0, 3),
      b_real ~ dnorm(0, 3),
      b_spur ~ dnorm(0, 3),
      sigma ~ dlnorm(0, 3)),
    data = Sim5.15Data
  ) 

m_real <-
  map(
    alist(
      y ~ dnorm(mu, sigma),
      mu <- a + b_real * x_real, 
      a ~ dnorm(0, 3),
      b_real ~ dnorm(0, 3),
      sigma ~ dlnorm(0, 3)),
    data = Sim5.15Data
  ) 

m_spur <- 
  map(
    alist(
      y ~ dnorm(mu, sigma),
      mu <- a + b_spur * x_spur,
      a ~ dnorm(0, 3),
      b_spur ~ dnorm(0, 3),
      sigma ~ dlnorm(0, 3)),
    data = Sim5.15Data
  )
lapply(list(both = m_both, real = m_real, spur = m_spur), precis)

## $both
##         Mean StdDev  5.5% 94.5%
## a      -0.29   0.10 -0.45 -0.13
## b_real  1.28   0.15  1.04  1.51
## b_spur -0.21   0.12 -0.40 -0.03
## sigma   0.98   0.07  0.87  1.09
## 
## $real
##         Mean StdDev  5.5% 94.5%
## a      -0.30   0.10 -0.46 -0.14
## b_real  1.06   0.10  0.91  1.22
## sigma   0.99   0.07  0.88  1.10
## 
## $spur
##         Mean StdDev  5.5% 94.5%
## a      -0.45   0.13 -0.66 -0.25
## b_spur  0.54   0.10  0.39  0.70
## sigma   1.28   0.09  1.14  1.42

gf_point(y ~ x_real, data = Sim5.15Data) %>%
  gf_coefline(coef = coef(m_real))

## Warning in gf_coefline(., coef = coef(m_real)): Ignoring all but first two
## values of coef.

gf_point(y ~ x_spur, data = Sim5.15Data) %>%
  gf_coefline(coef = coef(m_spur))

## Warning in gf_coefline(., coef = coef(m_spur)): Ignoring all but first two
## values of coef.

The milk data set

Source

Comparative primate milk composition data, from Table 2 of Hinde and Milligan. 2011. Evolutionary Anthropology 20:9-23.

Variables

• clade: Broad taxonomic group

• species: Species name

• kcal.per.g: Kilocalories per gram of milk

• perc.fat: Percent fat

• perc.protein: Percent protein

• perc.lactose: Percent lactose

• mass: Body mass of mother, in kilograms

• neocortex.perc: Percent of brain mass that is neocortex

R code 5.16

data(milk)

R code 5.17

# This fails.
m5.5 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + bn * neocortex.perc,
    a ~ dnorm(0, 100),
    bn ~ dnorm(0, 1),
    sigma ~ dunif(0, 1)
  ),
  data = milk
)

## Error in map(alist(kcal.per.g ~ dnorm(mu, sigma), mu <- a + bn * neocortex.perc, : initial value in 'vmmin' is not finite
## The start values for the parameters were invalid. This could be caused by missing values (NA) in the data or by start values outside the parameter constraints. If there are no NA values in the data, try using explicit start values.

R code 5.18

# Here's why the previous chunk failed -- missing data!
favstats( ~ neocortex.perc, data = milk)

##    min    Q1 median    Q3  max     mean       sd  n missing
##  55.16 64.54  68.85 71.26 76.3 67.57588 5.968612 17      12

R code 5.19

# Let's grab just the rows that have no missing data.
MilkCC <- milk %>% filter(complete.cases(.))
favstats(~ neocortex.perc, data = MilkCC)

##    min    Q1 median    Q3  max     mean       sd  n missing
##  55.16 64.54  68.85 71.26 76.3 67.57588 5.968612 17       0

R code 5.20

m5.5 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + bn * neocortex.perc,
    a ~ dnorm(0, 100),
    bn ~ dnorm(0, 1),
    sigma ~ dunif(0, 1)
  ),
  data = MilkCC
)

R code 5.21

precis(m5.5, digits = 3)

##        Mean StdDev   5.5% 94.5%
## a     0.353  0.471 -0.399 1.106
## bn    0.005  0.007 -0.007 0.016
## sigma 0.166  0.028  0.120 0.211

R code 5.22

coef(m5.5)["bn"] * (76 - 55)

##        bn 
## 0.0945674

R code 5.23

pred.data <- data.frame(neocortex.perc = 50:80)
mu <- link(m5.5, data = pred.data, n = 1e4)

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

m5.5.pred <-
  pred.data %>%
  mutate(
    mu.mean = apply(mu, 2, mean),
    mu.PIlo = apply(mu, 2, PI)[1,],
    mu.PIhi = apply(mu, 2, PI)[2,]
  )

gf_point(kcal.per.g ~ neocortex.perc + label:species, data = MilkCC,
         color = rangi2) %>%
  gf_line(mu.mean ~ neocortex.perc, data = m5.5.pred) %>%
  gf_ribbon(mu.PIlo + mu.PIhi ~ neocortex.perc, data = m5.5.pred,
            alpha = 0.2) %>%
plotly::ggplotly()

## Warning: Ignoring unknown aesthetics: label

### R code 5.24

MilkCC <- 
  MilkCC %>% 
  mutate(log.mass = log(mass))

R code 5.25

m5.6 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + bm * log.mass,
    a ~ dnorm(0, 100),
    bm ~ dnorm(0, 1),
    sigma ~ dunif(0, 1)
  ),
  data = MilkCC
)
precis(m5.6)

##        Mean StdDev  5.5% 94.5%
## a      0.71   0.05  0.63  0.78
## bm    -0.03   0.02 -0.06  0.00
## sigma  0.16   0.03  0.11  0.20

R code 5.26

m5.7 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + bn * neocortex.perc + bm * log.mass,
    a ~ dnorm(0, 100),
    bn ~ dnorm(0, 1),
    bm ~ dnorm(0, 1),
    sigma ~ dunif(0, 1)
  ),
  data = MilkCC,
  start = list(a = 0, bn = 0, bm = 0, sigma = 0.50)
)
precis(m5.7)

##        Mean StdDev  5.5% 94.5%
## a     -1.08   0.47 -1.83 -0.34
## bn     0.03   0.01  0.02  0.04
## bm    -0.10   0.02 -0.13 -0.06
## sigma  0.11   0.02  0.08  0.15

R code 5.27

m5.7.pred <-
  data_frame(
    neocortex.perc = 50:80,
    log.mass = mean(log(MilkCC$mass))
  )
m5.7.link <- link(m5.7, data = m5.7.pred, n = 1e4)

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

m5.7.pred <-
  m5.7.pred %>% 
  mutate(
    mu.mean = apply(mu, 2, mean),
    mu.PIlo = apply(mu, 2, PI)[1,],
    mu.PIhi = apply(mu, 2, PI)[2,]
  )

gf_point(kcal.per.g ~ neocortex.perc + label:species, data = MilkCC) %>%
gf_ribbon(mu.PIlo + mu.PIhi ~ neocortex.perc, data = m5.7.pred, alpha = 0.2) %>%
gf_line(mu.mean ~ neocortex.perc, data = m5.7.pred, alpha = 0.5) %>%
plotly::ggplotly()

## Warning: Ignoring unknown aesthetics: label

R code 5.28

# simulating a masking relationship
# n = number of cases
# rho = correlation between two variables x_pos and x_neg
sim_masking <- function(n = 100, rho = 0.7) {
  data_frame(
    x_pos = rnorm(n),
    x_neg = rnorm(n, rho * x_pos, sd = sqrt(1 - rho^2)),
    y = rnorm(n, x_pos - x_neg, sd = 1)   # y equally associated to each var
  )
}
MaskingData <- sim_masking()
splom(MaskingData)        # splom = scatter plot matrix (lattice)
GGally::ggpairs(MaskingData)   # ggplot2 version

R code 5.29

sim_legs <- function(n = 100) {
  data_frame(
    height = rnorm(n, 10, 2),    # units = ??
    leg_prop = runif(n, 0.4, 0.5),
    leg_left = leg_prop * height + rnorm(n, 0, 0.2),
    leg_right = leg_prop * height + rnorm(n, 0, 0.2)
  )
}
LegHeight <- sim_legs()

R code 5.30

m5.8 <- map(
  alist(
    height ~ dnorm(mu, sigma),
    mu <- a + bl * leg_left + br * leg_right,
    a ~ dnorm(10, 100),
    bl ~ dnorm(2, 10),
    br ~ dnorm(2, 10),
    sigma ~ dunif(0, 10)
  ),
  data = LegHeight
)
precis(m5.8)

##       Mean StdDev 5.5% 94.5%
## a     1.44   0.40 0.81  2.08
## bl    1.42   0.26 1.00  1.84
## br    0.49   0.26 0.06  0.91
## sigma 0.72   0.05 0.64  0.81

R code 5.31

GGally::ggpairs(LegHeight)
plot(precis(m5.8))

R code 5.32

m5.8.post <- extract.samples(m5.8)

gf_point(bl ~ br, data = m5.8.post, col = rangi2, alpha = 0.1, pch = 16)

R code 5.33

# clunky for now -- on the fly addition not working in gf functions (yet)
m5.8.post <-
  m5.8.post %>%
  mutate(sum_bl_br = bl + br)
gf_dens( ~ sum_bl_br, data = m5.8.post) +
  labs(x = "bl + br")

R code 5.34

m5.9 <- map(alist(
  height ~ dnorm(mu, sigma),
  mu <- a + bl * leg_left,
  a ~ dnorm(10, 100),
  bl ~ dnorm(2, 10),
  sigma ~ dunif(0, 10)
),
data = LegHeight)
precis(m5.9)

##       Mean StdDev 5.5% 94.5%
## a     1.55   0.40 0.91  2.18
## bl    1.88   0.08 1.75  2.02
## sigma 0.74   0.05 0.65  0.82

R code 5.35

library(rethinking)
data(milk)

R code 5.36

# kcal.per.g regressed on perc.fat
m5.10 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + bf * perc.fat,
    a ~ dnorm(0.6, 10),
    bf ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = milk
)

# kcal.per.g regressed on perc.lactose
m5.11 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + bl * perc.lactose,
    a ~ dnorm(0.6, 10),
    bl ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = milk
)

lapply(list(fat = m5.10, lactose = m5.11), precis, digits = 3)

## $fat
##        Mean StdDev  5.5% 94.5%
## a     0.301  0.036 0.244 0.358
## bf    0.010  0.001 0.008 0.012
## sigma 0.073  0.010 0.058 0.089
## 
## $lactose
##         Mean StdDev   5.5%  94.5%
## a      1.166  0.043  1.098  1.235
## bl    -0.011  0.001 -0.012 -0.009
## sigma  0.062  0.008  0.049  0.075

R code 5.37

m5.12 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + bf * perc.fat + bl * perc.lactose,
    a ~ dnorm(0.6, 10),
    bf ~ dnorm(0, 1),
    bl ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = milk
)
lapply(list(fat = m5.10, lactose = m5.11, `fat + lactose` = m5.12), precis, digits = 3)

## $fat
##        Mean StdDev  5.5% 94.5%
## a     0.301  0.036 0.244 0.358
## bf    0.010  0.001 0.008 0.012
## sigma 0.073  0.010 0.058 0.089
## 
## $lactose
##         Mean StdDev   5.5%  94.5%
## a      1.166  0.043  1.098  1.235
## bl    -0.011  0.001 -0.012 -0.009
## sigma  0.062  0.008  0.049  0.075
## 
## $`fat + lactose`
##         Mean StdDev   5.5%  94.5%
## a      1.007  0.200  0.688  1.327
## bf     0.002  0.002 -0.002  0.006
## bl    -0.009  0.002 -0.013 -0.005
## sigma  0.061  0.008  0.048  0.074

R code 5.38

GGally::ggpairs(milk %>% select(kcal.per.g, perc.fat, perc.lactose))

R code 5.39

cor(milk$perc.fat, milk$perc.lactose)

## [1] -0.9416373

# alternative syntax using mosaic package
cor(perc.fat ~ perc.lactose, data = milk)

## [1] -0.9416373

R code 5.40

library(rethinking)
data(milk)
collinearity_sim <- 
  function(rho = 0.9, data = milk) {
    data <- 
      data %>% 
      mutate(
        x = rnorm(nrow(data), mean = rho * perc.fat,
               sd = sqrt((1 - rho^2) * var(perc.fat)))
      )
    model = lm(kcal.per.g ~ perc.fat + x, data = data)
    sqrt(diag(vcov(model)))[2] # stddev of parameter
  } 
collinearity_sim <-
  Vectorize(collinearity_sim, "rho")

SimData5.40 <-
  expand.grid(
    r = seq(from = 0, to = 0.99, by = 0.01),
    rep = 1:100) %>%
  mutate(
    sd = collinearity_sim(r)
  ) 

gf_point(sd ~ r, data = SimData5.40, alpha = 0.01) %>%
  gf_spline(sd ~ r)

R code 5.41

# simulate data where growth is inhibited by fungus, which is inhibit by soil treatments
# number of plants
SimData5.13 <-
  expand.grid(
    treatment = c(0, 1),
    rep = 1:50              # 50 plants in each treatment group
  ) %>% 
  mutate(
    height0 = rnorm(100, 10, 2),  # initial heights of plants
    # fungus grows in half of control group and 10% of treatment group
    fungus = rbinom(100, size = 1, prob = 0.5 - 0.4 * treatment),
    # mean growth is 5 without fungus, 2 with fungus
    height1 = height0 + rnorm(100, 5 - 3 * fungus)
  )

R code 5.42

If we use treatment and fungus to predict growth, treatment appears to have no effect. But we know it does (since we simulated it that way). The impact of treatment is masked by the use of fungus, since the way treatment affects growth is by inhibiting fungus.

m5.13 <- map(
  alist(
    height1 ~ dnorm(mu, sigma),
    mu <- a + bh * height0 + bt * treatment + bf * fungus,
    a ~ dnorm(0, 100),
    c(bh, bt, bf) ~ dnorm(0, 10),
    sigma ~ dunif(0, 10)
  ),
  data = SimData5.13
)
precis(m5.13)

##        Mean StdDev  5.5% 94.5%
## a      5.13   0.54  4.26  6.00
## bh     1.01   0.05  0.93  1.09
## bt    -0.28   0.25 -0.67  0.12
## bf    -2.94   0.27 -3.37 -2.52
## sigma  1.06   0.07  0.94  1.18

R code 5.43

If we remove fungus, we can see the effect of treatment.

m5.14 <- map(
  alist(
    height1 ~ dnorm(mu, sigma),
    mu <- a + bh * height0 + bt * treatment,
    a ~ dnorm(0, 100),
    c(bh, bt) ~ dnorm(0, 10),
    sigma ~ dunif(0, 10)
  ),
  data = SimData5.13,
  start = list(a = 2, bh = 1, bt = 1, sigma = 1)
)
precis(m5.14)

##       Mean StdDev 5.5% 94.5%
## a     2.84   0.75 1.64  4.04
## bh    1.07   0.07 0.95  1.18
## bt    1.19   0.31 0.69  1.69
## sigma 1.57   0.11 1.39  1.75

R code 5.44

data(Howell1)

R code 5.45

m5.15 <- 
  map(
    alist(
      height ~ dnorm(mu, sigma),
      mu <- a + bm * male,
      a ~ dnorm(178, 100),
      bm ~ dnorm(0, 10),
      sigma ~ dunif(0, 50)
    ),
    data = Howell1)
precis(m5.15)

##         Mean StdDev   5.5%  94.5%
## a     134.83   1.59 132.29 137.38
## bm      7.28   2.28   3.63  10.93
## sigma  27.31   0.83  25.99  28.63

R code 5.46

m5.15.post <- 
  extract.samples(m5.15) %>% 
  mutate(
    mu.male = a + bm
  )
PI(m5.15.post$mu.male)

##       5%      94% 
## 139.4116 144.8732

R code 5.47

m5.15b <- 
  map(
    alist(
      height ~ dnorm(mu, sigma),
      mu <- af * (1 - male) + am * male,
      af ~ dnorm(178, 100),
      am ~ dnorm(178, 100),
      sigma ~ dunif(0, 50)
    ),
    data = Howell1,
    start = list(af = 178, am = 178, sigma = 3)
  )

R code 5.48

data(milk)
tally(~ clade, data = milk)

## clade
##              Ape New World Monkey Old World Monkey    Strepsirrhine 
##                9                9                6                5

R code 5.49 - 50

milk <-
  milk %>% 
  mutate(
    clade.NWM = ifelse(clade == "New World Monkey", 1, 0),
    clade.OWM = ifelse(clade == "Old World Monkey", 1, 0),
    clade.S = ifelse(clade == "Strepsirrhine", 1, 0)
  )

R code 5.51

m5.16 <- map(
  alist(
    kcal.per.g ~ dnorm(mu, sigma),
    mu <- a + b.NWM * clade.NWM + b.OWM * clade.OWM + b.S * clade.S,
    a ~ dnorm(0.6, 10),
    b.NWM ~ dnorm(0, 1),
    b.OWM ~ dnorm(0, 1),
    b.S ~ dnorm(0, 1),
    sigma ~ dunif(0, 10)
  ),
  data = milk
)
precis(m5.16)

##        Mean StdDev  5.5% 94.5%
## a      0.55   0.04  0.49  0.61
## b.NWM  0.17   0.05  0.08  0.25
## b.OWM  0.24   0.06  0.15  0.34
## b.S   -0.04   0.06 -0.14  0.06
## sigma  0.11   0.02  0.09  0.14

R code 5.52

# sample posterior
m5.16.post <- 
  extract.samples(m5.16) %>%
  mutate(
    # compute averages for each category
    mu.ape = a,
    mu.NWM = a + b.NWM,
    mu.OWM = a + b.OWM,
    mu.S   = a + b.S
  )

# summarize using precis 
# computes mean, sd and PI for each variable in data frame
precis(m5.16.post, digits = 3)

##          Mean StdDev  |0.89 0.89|
## a       0.546  0.038  0.486 0.608
## b.NWM   0.168  0.053  0.084 0.255
## b.OWM   0.242  0.060  0.147 0.338
## b.S    -0.039  0.064 -0.141 0.062
## sigma   0.114  0.015  0.092 0.140
## mu.ape  0.546  0.038  0.486 0.608
## mu.NWM  0.714  0.038  0.654 0.774
## mu.OWM  0.788  0.046  0.711 0.858
## mu.S    0.507  0.051  0.431 0.593

R code 5.53

quantile( ~ (mu.NWM - mu.OWM), data = m5.16.post, probs = c(0.025, 0.5, 0.975))

##        2.5%         50%       97.5% 
## -0.19060738 -0.07364717  0.04551348

R code 5.54

milk <-
  milk %>% 
  mutate(
    clade_id = coerce_index(clade)
  )

R code 5.55

m5.16_alt <- 
  map(
    alist(
      kcal.per.g ~ dnorm(mu, sigma),
      mu <- a[clade_id],
      a[clade_id] ~ dnorm(0.6, 10),
      sigma ~ dunif(0, 10)
    ),
    data = milk)
precis(m5.16_alt, depth = 2)

##       Mean StdDev 5.5% 94.5%
## a[1]  0.55   0.04 0.48  0.61
## a[2]  0.71   0.04 0.65  0.78
## a[3]  0.79   0.05 0.71  0.86
## a[4]  0.51   0.05 0.43  0.59
## sigma 0.11   0.02 0.09  0.14

R code 5.56

m5.17 <- lm(y ~ 1 + x, data = d)
m5.18 <- lm(y ~ 1 + x + z + w, data = d)

R code 5.57

m5.19a <- lm(y ~ 1 + x, data = d)
m5.19b <- lm(y ~ x, data = d)

R code 5.58

m5.20 <- lm(y ~ 0 + x, data = d)
m5.21 <- lm(y ~ x - 1, data = d)

R code 5.59

m5.22 <- lm(y ~ 1 + as.factor(season), data = d)

R code 5.60

m5.23 <- lm(y ~ 1 + x + x2 + x3, 
            data = d %>% mutate(x2 = x^2, x3 = x^3))

R code 5.61

m5.24 <- lm(y ~ 1 + x + I(x^2) + I(x^3), data = d)

R code 5.62

data(cars)
glimmer(dist ~ speed, data = cars)

## alist(
##     dist ~ dnorm( mu , sigma ),
##     mu <- Intercept +
##         b_speed*speed,
##     Intercept ~ dnorm(0,10),
##     b_speed ~ dnorm(0,10),
##     sigma ~ dcauchy(0,2)
## )